##
# Replication Files for Legislative Resources, Corruption, and Incumbency
# British Journal of Political Science
# Shane Martin, Charles T. McClean, and Kaare W. Strom
##

#################
# Load Packages #
#################

library(tidyverse)
library(stargazer)
library(lfe)
library(interflex)
library(AER)
library(xtable)

###################
# Read in Dataset #
###################

cy <- readRDS("incumbency.RDS")

#############
# savepdf() #
#############

savepdf <- function(file, width=16, height=10){ 
  fname <- paste("",file,".pdf",sep="")
  pdf(fname, width=width/2.54, height=height/2.54,
      pointsize=10)
  par(mgp=c(0,0.45,0), tcl=-0.4, mar=c(3.3,3.6,1.1,1.1)) # cex = 0.8, mar 0 was 1.1
}

###########################################################
# FIGURE 1: Incumbency Rate in 68 Democracies (2000-2018) #
###########################################################

c <- cy %>%
  group_by(Country) %>%
  summarize(median_inc = median(IncumbencyRate)) %>%
  arrange(median_inc)

cy$Country <- factor(cy$Country, levels = rev(c$Country))
rm(c)

ggplot(cy, aes(x = Country, y = IncumbencyRate)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(name = "Incumbency Rate",
                   breaks = c(0, 25, 50, 75, 100),
                   limits = c(0,100)) +
  scale_x_discrete(name = "") +
  theme_bw()

ggsave("1_Boxplot.pdf",
       plot = last_plot(),
       device = "pdf",
       width = 6,
       height = 9,
       units = "in")

######################################
# Exclude Countries With Term Limits #
######################################

cy2 <- cy %>%
  filter(!(Country %in% c("COSTA RICA", "ECUADOR", "MEXICO", "BOLIVIA", "PHILIPPINES")))

##############################################################
# TABLE 1: Legislative Resources, Corruption, and Incumbency #
##############################################################

m1 <- felm(IncumbencyRate ~ LegOrg + Corruption | 0 | 0 | Country, data=cy2)
m2 <- felm(IncumbencyRate ~ LegOrg + Corruption + Bicameral + Presidential + Federal + ICPV + LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m3 <- felm(IncumbencyRate ~ LegOrg * Corruption | 0 | 0 | Country, data=cy2)
m4 <- felm(IncumbencyRate ~ LegOrg * Corruption + Bicameral + Presidential + Federal + ICPV + LogGDP + GDPPct | 0 | 0 | Country, data=cy2)

stargazer(m1, m2, m3, m4, type = "text",
          dep.var.labels = c("Incumbency Rate"),
          omit.stat = c("f", "ser", "adj.rsq"),
          digits=3,
          column.sep.width = "0pt",
          covariate.labels = c("Legislative Resources", "Corruption",
                               "Bicameral", "Presidential", "Federal",
                               "ICPV",
                               "log(GDP per Capita)", "GDP Growth",
                               "Legislative Resources x Corruption"),
          title="Legislative Resources, Corruption, and Incumbency") 

rm(m1,m2,m3,m4)

########################################################################################
# FIGURE 2: Marginal Effect of Legislative Resources on Incumbency by Corruption Level #
########################################################################################

savepdf("2_Marginal")
out <- interflex(Y = "IncumbencyRate", D = "LegOrg", X = "Corruption", data=cy2, estimator="linear", vartype="delta",
                 main = "", Z=c("Presidential", "Federal", "Bicameral", "ICPV", "GDPPct", "LogGDP"), cl="Country",
                 xlab = "Moderator: Corruption", ylab="Marginal Effect of Legislative Resources on Incumbency Rate",
                 cex.axis=0.5, cex.lab=0.5, full.moderate = FALSE, theme.bw = TRUE, show.grid = FALSE, ylim=c(-0.6,0.9)) 
out$figure
dev.off()

rm(out)

##############################################################
# TABLE 2: Legislative Resources, Corruption, and Incumbency #
##############################################################

m1 <- ivreg(IncumbencyRate ~ LegOrg + Corruption | log(WordCount) + Corruption, data=cy2)
m2 <- ivreg(IncumbencyRate ~ LegOrg + Corruption + Bicameral + Presidential + Federal + ICPV + LogGDP +
              GDPPct | log(WordCount) + Corruption + Bicameral + Presidential + Federal + ICPV +
              LogGDP + GDPPct, data=cy2)
m3 <- ivreg(IncumbencyRate ~ LegOrg * Corruption | log(WordCount) * Corruption, data=cy2)
m4 <- ivreg(IncumbencyRate ~ LegOrg * Corruption + Bicameral + Presidential + Federal + ICPV + LogGDP +
              GDPPct | log(WordCount) * Corruption + Bicameral + Presidential + Federal + ICPV +
              LogGDP + GDPPct, data=cy2)

summary(m1, vcov=vcovHC(m1, type="HC1"))
summary(m2, vcov=vcovHC(m2, type="HC1"))
summary(m3, vcov=vcovHC(m3, type="HC1"))
summary(m4, vcov=vcovHC(m4, type="HC1"))

rm(m1,m2,m3,m4)

############
# APPENDIX #
############

################################
# TABLE A3: Summary Statistics #
################################

vars <- c("IncumbencyRate", "LegOrg", "CommNo", "Staff", "CommMult", "LogSalary", "ICPV",
          "Bicameral", "Presidential", "Federal", "GDPPct", "LogGDP", "Corruption")
results <- as.data.frame(matrix(NA, ncol=6, nrow=13))
colnames(results) <- c("Variable", "Mean", "SD", "Min", "Max", "N")
subset <- cy2[,vars]
results[,1] <- c("Incumbency Rate", "Legislative Resources Index", "Number of Committees", "Staff per Legislator", "Multiple Committees",
                 "log(Salary/GNI per Capita)", "ICPV", "Bicameral", "Presidential", "Federal", 
                 "GDP Growth", "log(GDP per Capita)", "Corruption")
for(i in 1:length(results$Variable)){
  results[i,2] <- mean(subset[,i], na.rm=T)
  results[i,3] <- round(sd(subset[,i], na.rm=T),3)
  results[i,4] <- min(subset[,i], na.rm=T)
  results[i,5] <- max(subset[,i], na.rm=T)
  results[i,6] <- length(subset[,i])
}
results$SD <- as.character(results$SD)
results$SD <- paste("(",  results$SD, sep="")
results$SD <- paste(results$SD, ")", sep="")
print(xtable(results, digits = 2), include.rownames=FALSE)

rm(results, i, vars, subset)

###############################################################################
# TABLE A4: Legislative Resources, Corruption, and Incumbency (Fixed Effects) #
###############################################################################

m1 <- felm(IncumbencyRate ~ LegOrg + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | Year | 0 | Country, data=cy2)
m2 <- felm(IncumbencyRate ~ Corruption | Year + Country | 0 | Country, data=cy2)
m3 <- felm(IncumbencyRate ~ LegOrg * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | Year | 0 | Country, data=cy2)

stargazer(m1, m2, m3, type = "text",
          dep.var.labels = c("Incumbency Rate"),
          omit.stat = c("f", "ser", "adj.rsq"),
          digits=3,
          column.sep.width = "0pt",
          covariate.labels = c("Legislative Resources", "Corruption",
                               "Bicameral", "Presidential", "Federal",
                               "ICPV", 
                               "log(GDP per Capita)", "GDP Growth",
                               "Legislative Resources x Corruption"),
          title="Legislative Resources, Corruption, and Incumbency (Fixed Effects)")

rm(m1,m2,m3)

########################################################################################
# TABLE A5: Legislative Resources, Corruption, and Incumbency (Country-Level Analysis) #
########################################################################################

avg <- cy2 %>%
  group_by(Country) %>%
  summarize(IncumbencyRate = median(IncumbencyRate),
            LegOrg = median(LegOrg),
            Corruption = median(Corruption),
            Bicameral = median(Bicameral),
            Presidential = median(Presidential),
            Federal = median(Federal),
            ICPV = median(ICPV),
            LogGDP = median(LogGDP),
            GDPPct = median(GDPPct))

m1 <- felm(IncumbencyRate ~ LegOrg + Corruption | 0 | 0 | 0, data=avg)
m2 <- felm(IncumbencyRate ~ LegOrg + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | 0, data=avg)
m3 <- felm(IncumbencyRate ~ LegOrg * Corruption | 0 | 0 | 0, data=avg)
m4 <- felm(IncumbencyRate ~ LegOrg * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | 0, data=avg)

stargazer(m1, m2, m3, m4, type = "text",
          dep.var.labels = c("Incumbency Rate"),
          omit.stat = c("f", "ser", "adj.rsq"),
          digits=3,
          column.sep.width = "0pt",
          covariate.labels = c("Legislative Resources", "Corruption",
                               "Bicameral", "Presidential", "Federal",
                               "ICPV", 
                               "log(GDP per Capita)", "GDP Growth",
                               "Legislative Resources x Corruption"),
          title="Legislative Resources, Corruption, and Incumbency (Average Values)")

rm(m1,m2,m3,m4,avg)

##################################################################################
# TABLE A6: Legislative Resources, Corruption, and Incumbency (Index Components) #
##################################################################################

m1 <- felm(IncumbencyRate ~ CommNo + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m2 <- felm(IncumbencyRate ~ CommNo * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m3 <- felm(IncumbencyRate ~ Staff + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m4 <- felm(IncumbencyRate ~ Staff * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m5 <- felm(IncumbencyRate ~ CommMult + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m6 <- felm(IncumbencyRate ~ CommMult * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m7 <- felm(IncumbencyRate ~ LogSalary + Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)
m8 <- felm(IncumbencyRate ~ LogSalary * Corruption + Bicameral + Presidential + Federal + ICPV +
             LogGDP + GDPPct | 0 | 0 | Country, data=cy2)

stargazer(m1, m2, m3, m4, m5, m6, m7, m8, type = "text",
          dep.var.labels = c("Incumbency Rate"),
          omit.stat = c("f", "ser", "adj.rsq"),
          digits=2,
          column.sep.width = "0pt",
          covariate.labels = c("No. of Committees", "Staff per Legislator", "Multiple Committees", "log(Salary/GNI per Capita)",
                               "Corruption",
                               "Bicameral", "Presidential", "Federal",
                               "ICPV", 
                               "log(GDP per Capita)", "GDP Growth",
                               "No. of Committees x Corruption", "Staff per Legislator x Corruption", "Multiple Committees x Corruption", "log(Salary/GNI per Capita) x Corruption"),
          title="Legislative Resources, Corruption, and Incumbency (Index Components)") 

rm(m1,m2,m3,m4,m5,m6,m7,m8)

############################################
# FIGURE A1: Marginal Effect of Corruption #
############################################

savepdf("A1_Corruption")
out <- interflex(Y = "IncumbencyRate", D = "Corruption", X = "LegOrg", data=cy2, estimator="linear", vartype="delta",
                 main = "", Z=c("Presidential", "Federal", "Bicameral", "ICPV", "GDPPct", "LogGDP"), cl="Country",
                 xlab = "Moderator: Legislative Resources", ylab="Marginal Effect of Corruption on Incumbency Rate",
                 cex.axis=0.5, cex.lab=0.5, full.moderate = FALSE, theme.bw = TRUE, show.grid = FALSE, ylim=c(-1.1,0.3)) 
out$figure
dev.off()

rm(out)

####################################
# FIGURE A2: Raw Data for Figure 2 #
####################################

savepdf("A2_Raw")
interflex(Y = "IncumbencyRate", D = "LegOrg", X = "Corruption", data=cy2, estimator="raw", vartype="delta",
                 main = "", Z=c("Presidential", "Federal", "Bicameral", "ICPV", "GDPPct", "LogGDP"), cl="Country",
                 Ylabel = "Incumbency Rate", Xlabel = "Corruption", Dlabel = "Legislative Resources",
                 cex.axis=0.5, cex.lab=0.5, full.moderate = FALSE, theme.bw = TRUE, show.grid = FALSE, ylim=c(-0.6,0.9), ncols=3) 
dev.off()

#############################################
# FIGURE A3: Binning Estimator for Figure 2 # 
#############################################

savepdf("A3_Binning")
out <- interflex(Y = "IncumbencyRate", D = "LegOrg", X = "Corruption", data=cy2, estimator="binning", vartype="delta", nbins = 3,
                 main = "", Z=c("Presidential", "Federal", "Bicameral", "ICPV", "GDPPct", "LogGDP"), cl="Country",
                 xlab = "Moderator: Corruption", ylab="Marginal Effect of Legislative Resources on Incumbency Rate",
                 cex.axis=0.5, cex.lab=0.5, full.moderate = FALSE, theme.bw = TRUE, show.grid = FALSE, ylim=c(-0.6,0.9), bin.labs = FALSE) 
out$figure
dev.off()

rm(out)

###############################################################
# FIGURE A4: Legislative Word Count and Legislative Resources # 
###############################################################

cons <- cy2 %>%
  select(Country, LegOrg, WordCount) %>%
  unique()

savepdf("A4_IV")
ggplot(cons, aes(x=log(WordCount), y=LegOrg)) +
  geom_point() +
  geom_smooth(method=lm, se=TRUE, color="black") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_text(size = 8),
        axis.text.y = element_text(size = 8),
        axis.text.x = element_text(size = 8)) +
  labs(x = "log(Legislative Word Count)", y = "Legislative Resources Index") +
  ylim(c(0,100))
dev.off()

rm(cons)
